Breaking adsorption-energy scaling limitations of electrocatalytic nitrate reduction on intermetallic CuPd nanocubes by machine-learned insights

The electrochemical nitrate reduction reaction (NO3RR) to ammonia is an essential step toward restoring the globally disrupted nitrogen cycle. In search of highly efficient electrocatalysts, tailoring catalytic sites with ligand and strain effects in random alloys is a common approach but remains limited due to the ubiquitous energy-scaling relations. With interpretable machine learning, we unravel a mechanism of breaking adsorption-energy scaling relations through the site-specific Pauli repulsion interactions of the metal d-states with adsorbate frontier orbitals. The non-scaling behavior can be realized on (100)-type sites of ordered B2 intermetallics, in which the orbital overlap between the hollow *N and subsurface metal atoms is significant while the bridge-bidentate *NO3 is not directly affected. Among those intermetallics predicted, we synthesize monodisperse ordered B2 CuPd nanocubes that demonstrate high performance for NO3RR to ammonia with a Faradaic efficiency of 92.5% at −0.5 VRHE and a yield rate of 6.25 mol h−1 g−1 at −0.6 VRHE. This study provides machine-learned design rules besides the d-band center metrics, paving the path toward data-driven discovery of catalytic materials beyond linear scaling limitations.

N itrate (NO 3 − ) is one of the most common water pollutants from a variety of sources, including agricultural runoff, industrial wastewater discharges, and animal manures 1 . To remediate NO 3 − contamination, selective catalytic reduction systems have been actively pursued with the aim to harmonize the global nitrogen cycle (N-cycle) 2,3 involving the interconversion of dinitrogen (N 2 ) and reactive nitrogen species, e.g., ammonia (NH 3 ), nitrogen oxides (NO x ), and NO 3 − . In this regard, the electrocatalytic NO 3 − reduction reaction (NO 3 RR) to NH 3 with renewable electricity offers a practical path for restoring the disrupted N-cycle and, more importantly, a sustainable alternative to the energy-intensive Haber-Bosch process that results in 1-2% of global carbon dioxide (CO 2 ) emissions 4 . Of particular interest is the development of high-performance catalysts for handling high pH NO 3 − concentrates due to less formation of toxic byproduct NO x and the growing concern of removing NO 3 − from alkaline nuclear wastes 5 .
Copper (Cu) has demonstrated promise for catalyzing NO 3 RR in alkaline media with reasonable Faradaic efficiencies (FE) to NH 3 1 , although high overpotentials are needed. Single-crystal experiments showed that NO 3 RR on Cu is structure sensitive with (100)-oriented surface sites more active toward NH 3 formation than the (111) counterparts 5 . A series of random alloy electrocatalysts (e.g., CuNi 6 , CuRh 7 , and PtRu 8 ) have been synthesized for NO 3 RR with a general tradeoff of the partial current density and FE, arguably due to the ubiquitous adsorption-energy scaling relations. Many strategies have been visioned to circumvent such energy-scaling limitations 9 on catalytic performance, such as tuning strain 10 and ligand 11 , designing bifunctional 12 or molecular single-site catalysts 13 , and imposing nanoscopic confinement 14 . However, the lack of theoretical underpinning of site reactivity makes it difficult to implement those strategies. Therefore, it is imperative to develop theory-guided principles that are transformative in material design for the development of advanced catalytic systems, particularly for finding high-performance NO 3 RR electrocatalysts toward NH 3 . Ordered intermetallic alloys, with atomically ordered structures and well-defined compositions, have attracted extensive attention as excellent electrocatalysts for oxygen reduction [15][16][17] , small molecules oxidation 18,19 , and CO 2 reduction 20 . Compared with random alloys, the structural ordering of intermetallic nanocrystals endows them with unique electronic properties and chemical stability 21 . Nevertheless, the direct solution-phase synthesis of ordered intermetallic nanocrystals remains challenging and their structural effect on surface reaction kinetics is largely unexploited.
In this work, we develop a mechanistic understanding of NO 3 RR on Cu surfaces in alkaline media with grand-canonical density functional theory (DFT) calculations. With the adsorption energies of bridge-bidentate *NO 3 and hollow *N as reactivity descriptors, the volcano plot very well captures the known activity trends of pure metals. The Bayesian theory of chemisorption (Bayeschem) 22 , as an interpretable machine learning (ML) approach, then unravels the origin of linear scaling between *NO 3 and *N adsorption energies on metal surfaces and identified an intriguing mechanism to break the scaling by leveraging the site-specific Pauli repulsion of metal d-states with adsorbate frontier orbitals. The machine-learned insights point to the peculiar properties of (100)-oriented surface sites at ordered intermetallics of d metals with a body-centered cubic (bcc) structure (B2) in which a shortened interlayer spacing along the bcc {100} direction results in significant overlap of the subsurface metal d-orbitals with the hollow *N p-orbitals while the bridgebidentate *NO 3 is not directly influenced. We then synthesize ordered intermetallic B2 CuPd nanocubes terminated with (100) facets using a colloidal method. The CuPd catalyst exhibits superior NO 3 RR performance compared with Cu and Pd nanocubes in alkaline media, validating theoretical predictions. Specifically, these CuPd nanocubes demonstrate a NH 3 FE of 92.5% at −0.5 V vs. reversible hydrogen electrode (RHE) and a yield rate of 6.25 mol h −1 g −1 at −0.6 V vs. RHE in NO 3 RR. Furthermore, these B2 CuPd nanocubes demonstrate high stability over 12 h electrolysis in 1 M KNO 3 + 1 M KOH. Bayeschem models suggest that while the upshifted d-band center of Cu sites at CuPd nanocubes favors the bridge-bidentate *NO 3 adsorption, the hollow *N is destabilized due to a dominant role of Pauli repulsion from the subsurface Pd d-orbitals, facilitating the protonation of N-bonded species toward NH 3 . This study demonstrates the concept of combining interpretable ML with precision synthesis for designing catalytic systems that break adsorptionenergy scaling relations and circumvent the corresponding limitations on attainable catalytic performance.

Results
Structure-activity relationships of NO 3 RR on metal surfaces. DFT calculations using the Vienna Ab initio Simulation Package (VASP) were performed to probe reaction pathways of NO 3 RR on metal surfaces. In the NO 3 RR literature, there are several suggested reaction pathways toward NH 3 formation depending on the electrolyte pH and catalysts 5,23,24 , and there is no consensus regarding the critical intermediates and governing factors. Supplementary Figure 1 shows the free energy diagram of different reaction pathways on Cu(100) and Cu(111) at 0 V vs. RHE in alkaline conditions. We chose the pathway *NO 3 → *NO 2 → * NO → *NHO → *NH 2 O → *NH 2 OH → *NH 2 → NH 3 because it is the most thermodynamically favorable one on both Cu(100) and Cu(111). Figure 1a, b shows the structures of reaction intermediates and free energy profiles of NO 3 RR to NH 3 on Cu(100) and Cu(111) at 0 V vs. RHE calculated from the grand-canonical DFT approach (see Supplementary Tables 1-3 for free energy corrections and source data). In alkaline media, the NO 3 RR is assumed to follow a series of deoxidation steps to form *NO, and its further reduction to NH 3 which likely goes through *NHO to adsorbed hydroxylamine (*NH 2 OH) 5 . It is generally accepted that NH 2 OH can only be transiently observed 25 , which is then readily reduced to NH 3 . At 0 V vs. RHE, the overpotential (E 0 = +0.69 V vs. RHE at pH 14) drives the removal of N-bonded species through electrochemical steps. However, the adsorption of negatively charged *NO 3 species (~−1.0 e for both surfaces from Bader charge analysis) is thermodynamically uphill, which is consistent with the observation that the adsorption of NO 3 − ions is rate-limiting for NO 3 RR on Cu at highly reductive potentials. The more favorable adsorption of NO 3 − on Cu(100) than that on Cu(111) explains an earlier onset potential for NO 3 − reduction on Cu(100) 5 . Moreover, *NO 2 formation from NO 2 − is more endergonic on Cu(111) than that on Cu(100), resulting in a slower re-adsorption of NO 2 − and its further reduction on (111)-type sites as observed 24 .
To understand the activity trends of NO 3 RR across elemental metals, we have calculated the free formation energies of reaction intermediates on the (100) and (111) facets of late transition and noble metals at the most stable configurations. Linear adsorptionenergy scaling relations of reaction intermediates with *NO 3 and *N as descriptor species ( Supplementary Fig. 2), were used to develop the activity volcano plot (Fig. 1c). The adsorption energies of *N at the four-fold hollow site and *NO 3 at the bridge site were used because of the simplicity and accuracy for capturing the thermodynamic stability of reaction intermediates. We used the highest free energy change of all reaction steps at 0 V vs. RHE (close to experimentally measured onset potentials) to characterize the theoretical activity of (100)-and (111)-like metal sites. Although many promising approaches have been proposed to calculate electrochemical barriers 26-28 , very few benchmarks are available. The choice of thermodynamics-based descriptors instead of full kinetic analysis with explicitly computed activation barriers is deliberate for capturing general activity trends and guiding experimental design. The predicted activity trend of Cu(100) > Cu(111) is consistent with experimental measurements 5 . The volcano map also suggests that Cu is more active than Pd, Ag, and Au because of a favorable *NO 3 formation and facile removal of N-bonded species. It is notable that there is a strong linear scaling of adsorption energies between *NO 3 and *N for late transition metals with a slope of 1.19 (R 2 : 0.88), close to the theoretical slope of 1.5 29 . Coinage metals (Cu, Ag, and Au) typically follow energy scaling with a different slope due to repulsive interactions from the fully occupied metal d-band 30 . Nevertheless, the desire is to break such scaling relations toward the top of the activity volcano, i.e., finding the optimal catalysts that adsorb *NO 3 stronger and *N weaker than Cu(100). Although various strategies of harnessing internal and external factors in complex materials were envisioned 9,31 , there is no theoretical guidance for going beyond the linear adsorption-energy scaling relations that largely limit the attainable catalytic performance.
Physical insights from interpretable ML. To pinpoint the physical origin of the linear scaling relations between *NO 3 and *N adsorption energies on metal surfaces beyond the valency argument, we have employed a recently-developed Bayesian theory of chemisorption (Bayeschem) as an interpretable ML approach 15 . Built upon the d-band theory of chemisorption and Bayesian optimization by learning from ab initio adsorption properties, Bayeschem has been used for qualitatively understanding the nature of chemical bonding and underlying electronic factors governing the trend of surface reactivity 15 . Figure 2a Figure 2b shows DFT-calculated and model-predicted density of states projected onto adsorbate frontier orbitals, with Cu(100) as an example. For *N, both p xy and p z orbitals contribute to the adsorption energy with clearly captured bonding and antibonding states. In comparison, the HOMO of *NO 3 has antibonding states pinning across the Fermi level, while the LUMO is too high in energy to be occupied, thus forming a Lorenzian-shaped resonance state (Fig. 2b). By varying electronic factors in the Newns-Anderson model Hamiltonian 15 (using Cu as a reference), we are showing the Bayschem-predicted change of *NO 3 and *N adsorption energies in Figs. 2c, d for (100)-and (111)-like surface sites, respectively. As the d-band center of surface sites shifts around in response to a perturbation, for example, by +1.5 to −1.5 eV, the adsorption energies of *NO 3 and *N change side by side, i.e., a higher (lower) d-band center leads to a more (less) favorable interaction, irrespective of surface termination (Fig. 2c, d). Interestingly, with the increase of the interatomic coupling strength V 2 ad , *N adsorption at the (100) hollow site becomes stronger first due to the depopulation of adsorbate-metal antibonding states but weakens as the change of V 2 ad is larger than 4.0 eV 2 where Pauli repulsion becomes dominant. This nonlinear correlation results from a complex interplay of the orbital hybridization and orthogonalization, both of which depend on the coupling strength V 2 ad , albeit in different slopes  ( Supplementary Fig. 7a). An increase of the V 2 ad for Cu(111) leads to a monotonically strengthened adsorption of *N until its value becomes unrealistically large ( Supplementary Fig. 7b). Structural analysis of 8 site motifs (Supplementary Table 4 and Fig. 8) showed that the increase of interatomic coupling can be more drastic for the B2 (bcc) intermetallic structures because the subsurface metal-ligand is close to the hollow *N, as illustrated in Fig. 2e. The Bayeschem model suggests that the (100)-facet of B2 structural motifs can potentially break adsorption-energy scaling relations. With the theoretical guidance from the Bayeschem model, we show the DFT-calculated adsorption energies of *NO 3 and *N on (100)-terminated B2 intermetallics and on randomly sampled (100)-terminated face-centered cubic (fcc) intermetallics from the Materials Project (Fig. 2f). The randomly sampled (100)-surfaces follow linear scaling relations as suggested by the Bayeschem. Compared with Cu(100), (100)-terminated B2 CuPd, ZnRh, and ZnCu structures (the first element denotes the surface metal) are predicted to exhibit reactivity properties beyond the scaling relations, i.e., a stronger *NO 3 adsorption and a weaker *N adsorption (Fig. 2g). Zn-terminated ones are likely not stable in NO 3 RR operating conditions, and thus not considered for further studies. Cu-terminated CuPd (100) and (111) surfaces have been previously shown to be more stable than the Pd-termination using DFT-calculated surface energies 32,33 . DFTcalculated free energy diagrams of the full reaction pathway predict that the activity metric ΔG max at CuPd(100) is~0.10 eV  lower than that at Cu(100), with both surfaces limited by *NO 3 adsorption ( Supplementary Fig. 9). Although DFT calculations at the GGA-PBE level have an often-quoted error of ±0.2 eV for adsorption energies at metal surfaces, the relative error across similar systems is expected to be much smaller and the qualitative prediction is reliable in the methodology. To validate the physical understanding of chemical bonding attained from machine learning, we have performed a detailed electronic structure analysis of Cu(100) and CuPd(100) in Supplementary Fig. 10. As suggested by the Bayeschem ML models, the higher d-center of site Cu atoms at CuPd stabilizes *NO 3 while the larger interatomic coupling from the subsurface Pd ligand destabilizes the *N, realizing an independent tuning of the *N and *NO 3 binding energies to a certain extent. Guided by the physical insights obtained from interpretable ML instead of explicitly exploring the intermetallic design space (requiring at least 7226 DFT calculations, see Supplementary Table 4), we greatly speed up the screening process and quickly narrow down the B2 CuPd intermetallic as the target for synthesis.
Synthesis and structural characterizations of monodisperse intermetallic CuPd nanocubes. High-temperature annealing is usually involved in the synthesis of ordered intermetallic alloys in order to promote metal atom rearrangement and d-d orbital hybridization 18,[26][27][28] . Such a high-temperature process involves solid-state annealing and results in the sintering and aggregation of nanocrystals. Solution-phase synthesis of ordered intermetallic nanoparticles is more desirable and has been reported in the PdCu 34 and AuCu 21,35 systems. Unfortunately, it usually requires an additional step either electrochemically or through a seedmediated growth and diffusion process, and the produced intermetallic alloys are partially ordered [34][35][36] . 3a) and transmission electron microscopy (TEM) bright-field (BF) image (Fig. 3b) of the as-synthesized sample reveal a uniform cubic morphology with a size of 50 ± 4 nm. The TEM BF image of intermetallic CuPd nanocubes shows a brighter contrast at the center and a darker contrast at the edge, suggesting atom enrichment at the surface or corners. The atomic ordering in the CuPd nanocube can be directly visualized using Z contrast in the atomicresolution high-angle annular dark-field scanning transmission electron microscopy (HAADF-STEM) images (Fig. 3c, d). The corresponding Fourier-transform (FT) pattern of the CuPd nanocube in Fig. 3c inset and the atomic model overlay in Fig. 3d indicated that the ordered CuPd nanocube is a single crystalline with a B2 intermetallic cubic structure. A HAADF-STEM image (Fig. 3c) and a high-resolution TEM (HRTEM) image (Supplementary Fig. 16a) show that the particles are single-crystalline with a lattice spacing of 2.95 Å, which can be assigned to the (100) planes of the ordered CuPd intermetallic phase. The selected area electron diffraction (SAED) pattern of the CuPd nanocubes depicted in Supplementary Fig. 16b shows bright concentric rings that can be assigned to the (100), (110), and (200) planes of the ordered CuPd intermetallic phase, respectively. Moreover, as shown in Supplementary Fig. 17, the alternating intensity profile in the corresponding HAADF line profile further confirms the Cu/ Pd atomic ordering within the CuPd nanocube. To get further insight into the distribution of Pd and Cu in the as-synthesized CuPd nanocubes, elemental analysis was carried out. Figure 3e shows a HAADF-STEM image and the corresponding elemental maps of a representative CuPd nanocube. We observed that Cu and Pd are homogeneously distributed across the CuPd nanocube, which is consistent with a B2 intermetallic structure. As shown in Fig. 3f, the powder X-ray diffraction (XRD) pattern of the ordered CuPd nanocube demonstrates that the distinct peaks are completely consistent with the ordered CuPd B2 intermetallic phase (ICDD No. 01-078-4406) 17 . The appearance of characteristic superlattice peaks at 2θ = 30°and 43°c onfirms the structure of the B2 intermetallic phase. The composition of the ordered CuPd nanocubes was further investigated by using energy-dispersive X-ray spectroscopy (EDX) and inductively coupled plasma mass spectrometry (ICP-MS), and the obtained consistent results suggested that the Pd/Cu molar ratio is about 1:1. Only Pd and Cu could be detected in the EDX spectrum except for the carbon signal which comes from the carbon-coated TEM grid ( Supplementary Fig. 18). The electronic interactions between Pd and Cu were comprehensively investigated via multiple characterization techniques. X-ray photoelectron spectroscopy (XPS) measurements show that the binding energy of Pd 3d core levels upshifts by~0.68 eV versus pure Pd nanocubes (Fig. 4a), while the binding energy of Cu 2p core levels decreases by~0.80 eV versus pure Cu nanocubes (Fig. 4b), indicating the existence of charge transfer between Pd and Cu 37,38 . This charge transfer is further confirmed by the shift of absorption edge to the higher energy direction in X-ray absorption near-edge spectroscopy (XANES) of the Pd K-edge ( Fig. 4c and inset) and the Cu K-edge ( Supplementary Fig. 19a). Figure 4d and Supplementary Fig. 20 present Fouriertransformed Pd K-edge extended X-ray absorption fine structure (EXAFS) spectra of the ordered CuPd nanocubes as well as the reference (Pd foil). As shown in Fig. 4d, in comparison with Pd foil, the ordered CuPd nanocubes exhibit shorter interatomic distance R Pd-Cu(Pd) than that of cubic close-packed Pd, providing a structural basis for employing the Pauli repulsion to weaken the *N adsorption to promote NO 3 RR catalysis. Supplementary  Fig. 19b presents Fourier-transformed Cu K-edge EXAFS spectra of the ordered CuPd nanocubes and the reference (Cu foil). In comparison with Cu foil, the ordered CuPd nanocubes exhibit slightly longer interatomic distance R Cu-Pd(Cu) than that of cubic close-packed Cu. The shift of the first nearest coordination peaks of the ordered CuPd nanocubes demonstrates the slight change in interatomic distances (Fig. 4d). The fitting results of ordered CuPd nanocubes and Pd foil are listed in Supplementary Table 5, showing the smaller coordinate number of Pd atoms in ordered B2 CuPd nanocubes than that in bulk Pd foil, which is consistent with the theoretical coordination number of 8 for B2 structure and 12 for fcc structure. Moreover, the specific PdCu(Pd) bond length decreases from 2.74 Å for Pd foil to 2.62 Å for ordered CuPd nanocubes. Because of the reduced bond distance of CuPd, the B2-ordered intermetallic CuPd nanocubes provide a platform for leveraging Pauling repulsion to go beyond the energy-scaling relations discussed in our computational results.
In order to understand the formation mechanism of the ordered CuPd nanocubes, the nucleation and growth process were monitored through time-dependent experiments. The aliquots taken from the synthesis process at different time intervals were characterized by TEM to capture the morphological evolution of the ordered CuPd nanocubes (Supplementary Fig. 21). In the beginning, when the reaction temperature was raised to 250°C for 1 min, only small spherical nanoparticles and some aggregates were observed. As the reaction time prolonged, nanocubes emerged. By the reaction time of 20 min, CuPd nanocubes were formed. The ICP results of the samples obtained at different reaction times suggest that Cu nanoparticles formed at the initial stage of the reaction and served as seeds for Pd to nucleate and grow (Supplementary Table 6). With the extension of reaction time, they grew into Cu-rich CuPd nanocubes. When the reaction was prolonged to 30 min, the ordered CuPd nanocubes formed. The final ordered CuPd structure with atom enrichment at the surface and corners could be attributed to the combination of galvanic replacement and Kirkendall effect 39,40 . This observation indicates that a multistage nucleation and growth process could promote the formation of ordered intermetallic structures, which was previously facilitated by tedious additional treatments (e.g., thermal annealing 18,41-43 , electrochemical methods 34 ). NO 3 RR on ordered intermetallic CuPd nanocubes. In order to evaluate the NO 3 RR performance of the ordered CuPd nanocubes, these nanocubes were loaded onto carbon black (Vulcan XC-72R) and treated in acetic acid to remove the surfactants according to previously reported methods 44,45 (Supplementary Fig. 22). As a comparison, Cu nanocubes and Pd nanocubes were synthesized. As shown in Supplementary Figs. 23, 24, the size of Cu nanocubes is about 35 nm, and the size of Pd nanocubes is about 12 nm. The XRD patterns show that both of them are fcc structures (Supplementary Figs. 23b, 24b). Electrochemical measurements were performed using a three-electrode system in a gas-tight H-cell separated by an ion-exchange membrane (Nafion 117) at room temperature. Pt foil and Ag/AgCl (3.5 M KCl) were used as the counter and reference electrodes, respectively. The electrocatalysts were deposited onto 1 cm 2 carbon fiber paper, leading to a metal loading of~0.2 mg cm −2 . The polarization curves were obtained by sweeping the potential from −0.5 to −1.7 V vs. Ag/AgCl at room temperature with a sweep rate of 20 mV s −1 in the Ar-saturated 1 M KOH + 1 M KNO 3 solution at room temperature. As shown in    Figure 5b shows that the partial current densities of NH 3 on ordered CuPd nanocubes are much higher than that on Pd nanocubes and Cu nanocubes, suggesting that the ordered B2 CuPd nanocubes are much more active than Pd nanocubes and Cu nanocubes for the NO 3 RR toward NH 3 . As shown in Supplementary Fig. 25, the ordered CuPd nanocubes show higher hydrogen evolution reaction (HER) activity than Cu nanocubes while Pd nanocubes perform best among all three tested. Cyclic voltammogram (CV) measurements at various scan rates (20,40, and 60 mV s −1 , etc.) were conducted in static solution to estimate the double-layer capacitance (C dl ) by sweeping the potential across the non-faradaic region −0.1-0 V vs. Ag/AgCl (Supplementary Fig. 26). The C dl for the ordered CuPd nanocubes is calculated to be 3.59 mF cm −2 , which is smaller than that of Cu nanocubes (4.93 mF cm −2 ) and Pd nanocubes (8.61 mF cm −2 ), due to the larger size of the ordered CuPd nanocubes. The intrinsic activities for NO 3 RR on ordered CuPd, Cu, and Pd nanocubes are evaluated by normalizing catalytic currents to ECSAs (Supplementary Fig. 27). The ECSA-normalized partial current densities of NH 3 on ordered CuPd nanocubes are much higher than that on Cu nanocubes and Pd nanocubes, indicating the intrinsic activity for NO 3 RR on ordered CuPd nanocubes is superior to those on Cu nanocubes and Pd nanocubes, consistent with DFT-predicted activity trends in Supplementary Fig. 9. Chronoamperometry (CA) measurements of catalysts were conducted at different potentials for 1 h in 1 M KOH + 1 M KNO 3 solution with continuous Ar bubbling at a rate of 20 standard cubic centimeters per minute (sccm) ( Supplementary  Fig. 28). The gas product was quantified by gas chromatography and only H 2 was identified from the competing HER. The colorimetric method using Nessler's reagent ( Supplementary   Fig. 29) was employed to detect the quantity of produced NH 3 and ion chromatography was employed to detect the quantity of produced NO 2 − (Supplementary Fig. 30). To avoid the loss of products from continuous Ar flow, a glass vial filled with 0.1 M HCl was set at the end of the outlet tube as the trap. The total NH 3 production yield was the summation of NH 3 in the electrolyte and 0.1 M HCl. The FEs and NH 3 yield rates of the electrocatalysts are shown in Fig. 5c, d. The ordered CuPd nanocubes demonstrated high selectivity toward NH 3 production from NO 3 RR with a FE of 92.5% at −0.5 V vs. RHE and a high yield rate of 6.25 mol h −1 g −1 at −0.6 V vs. RHE, outperforming most of the reported catalysts (Supplementary Table 7) 6,46-48 . The main byproduct of NO 3 RR on ordered CuPd nanocubes is NO 2 − , as detected and quantified by ion chromatography. The FE of NO 2 − starts from as high as 20.7% at −0.2 V vs. RHE, followed by a significant decrease to a minimal of~2.67% at −0.6 V vs. RHE. This suggested that NO 2 − could be an intermediate product and can be further reduced to NH 3 at a more negative potential, which is consistent with our theoretical results. Furthermore, in situ attenuated total reflectance surface-enhanced infrared absorption spectroscopy (ATR-SEIRAS) measurements were carried out to identify the intermediates of the NO 3 RR process. Supplementary Figure 31 shows the ATR-SEIRAS spectra results collected from ordered CuPd nanocubes during a CV cycle between 0.965 V and −0.835 V vs. RHE at 5 mV/s in 0.1 M KOH solution (Supplementary Fig. 31a) and 0.1 M KOH + 1 M KNO 3 solution (Supplementary Fig. 31b). According to previous reports, the absorption at 1645 cm −1 can be attributed to the H-O-H bending of water molecules 49 . The absorption at around 1365 cm −1 in Supplementary Fig. 31b is due to the adsorption of nitrate ions 50 . The intensity of the peak increased as the potential became more negative. There is no N=N stretching band at~2010 cm −1 appeared, which indicated that N 2 H x is not a reaction intermediate of the nitrate reduction on ordered CuPd nanocubes 51 .  To confirm the NO 3 RR performance, the control experiment was performed at −0.6 V vs. RHE for 1 h in 1 M KOH solution without KNO 3 . As shown in Fig. 5e, there is almost no NH 3 detected in the electrolyte. To confirm the produced NH 3 originated from the feeding nitrate solution, 15 N isotope labeling experiments were conducted. After electrolysis at −0.5 V vs. RHE for 1 h in 1 M K 15 NO 3 , no characteristic triple coupling peaks of 14 NH 4 + could be detected in the 1 H nuclear magnetic resonance ( 1 H NMR) spectra of the electrolyte, while only doublet peaks representing 15 NH 4 + were observed (Fig. 5f), indicating that the produced NH 3 entirely comes from the electroreduction of nitrate. To evaluate the long-term stability for future practical applications, we performed the durability test by running the chronoamperometry on the ordered CuPd nanocubes at −0.5 V vs. RHE for 12 h. As shown in Fig. 5g, the current density exhibits no appreciable decrease over 12 h of continuous operation, and the overall 12h FE of NH 3 is~85.1%. There is a slight decrease in the NO 3 RR electrocatalytic performance, possibly due to a small fraction of the catalyst falling off from the carbon paper. After the long-term electrolysis, the catalyst was removed from carbon paper and characterized by TEM and XPS. As shown in Supplementary Fig. 32, the cubic morphology is maintained, and there is no apparent aggregation observed. The HRTEM image (Supplementary Fig. 32c) shows clear lattice fringes with a lattice spacing of 2.95 Å, which can be assigned to the (100) planes of the ordered CuPd intermetallic phase. Moreover, the ICP-MS result shows that the Pd/Cu molar ratio is still close to 1:1. We also performed XRD measurement after the stability test ( Supplementary Fig. 33), which showed that the CuPd nanocubes still maintain the ordered B2 phase after long-term electrolysis. Furthermore, XPS analysis reveals that no obvious chemical state changes after the stability test ( Supplementary Fig. 34). These results illustrate the exceptional chemical and structural stability of the ordered CuPd nanocubes. Besides, electrochemical impedance spectroscopy (EIS) was recorded on these electrocatalysts to provide further insight into electrode kinetics ( Supplementary Fig. 35). Representative Nyquist plots show that the ordered CuPd nanocubes have a smaller charge transfer resistance (R ct ) than that of Cu nanocubes, indicating the fast Faradaic process and thus superior NO 3 RR kinetics for the ordered CuPd nanocubes, which stresses the importance of Pd in the intermetallic structure that modifies the electric conductivity of the ordered CuPd nanocubes.
Besides, in order to further validate the theoretical prediction of reactivity trends on pure metals in Fig. 1c, Cu spherical nanoparticles and Au particles (primary nanocubes) were synthesized (Supplementary Fig. 36). As shown in Supplementary  Fig. 37, Cu nanocubes with (100) facets show higher partial NH 3 Product from K NO 14 3 Fig. 5 Electrocatalytic NO 3 RR performance of the ordered CuPd nanocubes. a Linear scan voltammetry curves of ordered CuPd nanocubes, Cu nanocubes, and Pd nanocubes normalized to the geometric area. b Partial NH 3 current densities normalized to the geometric area. c FE of NH 3 at different potentials, The error bars correspond to the standard deviation from three independent measurements. d NH 3 yield rate of ordered CuPd nanocubes, Cu nanocubes, and Pd nanocubes at various potentials. e The yield rate of NH 3 with or without nitrate at −0.6 V vs. RHE. f NMR spectrum of the products generated during the electrocatalytic NO 3 RR with ordered CuPd nanocubes in 1 M K 15 NO 3 or 1 M K 14 NO 3 at −0.5 V vs. RHE. g Stability test by running the CA measurement on the ordered CuPd nanocubes at −0.5 V vs. RHE for 12 h. current density, FE and NH 3 yield rate than those of Cu nanoparticles, which indicates that Cu(100) has higher electrocatalytic activity for NO 3 RR to NH 3 than Cu(111) as DFTpredicted. The activities of NO 3 RR to NH 3 on nanocubes with (100) facets follow the trend Cu > Pd > Au, which is consistent with DFT-predicted activity trends in Fig. 1c.

Discussion
With electrocatalytic NO 3 − reduction to NH 3 on metal nanocatalysts as an example, we have demonstrated that interpretable ML of ab initio adsorption properties provides physical insights into the nature of chemical bonding that can be leveraged to break linear adsorption-energy scaling limitations of catalytic performance. Bayeschem models of two reactivity descriptors, i.e., bridge-bidentate *NO 3 and hollow *N, suggest that both adsorbates behave similarly in orbital hybridization upon a perturbation of the local electronic structure, e.g., d-band center. However, *N exhibits a more prevalent repulsion contribution on (100)-like sites than (111) as increasing interatomic coupling strengths, while the bridge-bidentate *NO 3 is not directly influenced due to large distances. These machine-learned insights can be leveraged for breaking linear scaling relationships, specifically at (100)terminated surface sites of B2 intermetallics, in which the layer separation is small compared to other site motifs of ordered intermetallics. A handful of B2 systems are predicted to be close to the top of the activity volcano plot because of a weakened *N binding and enhanced *NO 3 adsorption than Cu(100). Among those predicted, our synthesis strategy has enabled the direct, solution-phase formation of ordered intermetallic CuPd nanocubes that demonstrated highly efficient NO 3 RR. This study showcases a strategy for breaking the linear scaling relations on ordered intermetallic catalysts by harnessing the Pauli repulsion of the metal d-states with adsorbate frontier orbitals. Moreover, it highlights the benefit of interpretable ML and DFT calculations together with the synthesis of structurally controlled well-defined nanocrystals in suggesting and verifying the governing physical insights, providing a methodological basis for fine-tuning of electrocatalysts for improved efficiency.
Synthesis of ordered intermetallic CuPd nanocubes. In a modified procedure 52  Preparation of carbon-supported catalysts (20% loading). To prepare carbonsupported catalysts, the catalysts were loaded onto carbon black (Vulcan XC-72R) according to previously reported methods 44,45 . The hexane dispersion of 10 mg of intermetallic CuPd nanocubes was mixed with 40 mg carbon black and sonicated for 2 h. The product was collected by centrifugation. Afterwards, the catalysts were immersed in a mixture of 10 mL ethanol and 10 mL acetic acid for 10 h at 70°C to remove organic ligands on the surface of these CuPd nanocubes. The catalysts were washed three times with excessive ethanol and dried for 8 h in a vacuum oven at 60°C. The control samples of Cu and Pd nanocubes were loaded onto carbon via a similar approach.
Characterizations. XRD was performed on a Philips X' Pert PRO SUPER with Cu Kα (λ = 1.54056 Å). XPS was performed on a PHI Versa Probe III scanning XPS microscope using a monochromatic Al K-alpha X-ray source (1486.6 eV). The sample's morphology was characterized by SEM (Zeiss Supra 40) and TEM (EM-420). HRTEM, HAADF-STEM, X-EDS, and EDS mapping were conducted on a JEOL ARM 200CF equipped with an Oxford Instrument X-ray Energy Dispersive Spectrometer. The element contents of the products were determined by ICP-OES on a SPECTRO GENESIS ICP spectrometer. The gas product was quantified by gas chromatography (Agilent 7890B). The colorimetric method with Nessler's reagent on a UV-vis spectrophotometer (Agilent 3500) and ion chromatography (Metrohm Eco IC) was used to quantify the produced ammonia. Ion chromatography was also used to quantify the produced nitrite. The 1 H NMR signal was recorded on a Bruker 400 MHz system. The X-ray absorption spectra of Pd and Cu K-edges were obtained at the beamline 12-BM-B station of the Advanced Photon Source at Argonne National Laboratory. Both Pd and Cu K-edge XANES and EXAFS were measured under fluorescence mode by a Vortex ME4 detector. All XAS data analyses were performed with the Athena software package to extract XANES and EXAFS. Fourier-transform infrared spectroscopy (FTIR) was performed on an Agilent Cary 630.
Electrochemical measurements. Electrochemical measurements were carried out on a BioLogic electrochemical workstation. All measurements were performed in a gas-tight H-cell using a three-electrode system and an ion-exchange membrane (Nafion 117) at room temperature. Before testing, the Nafion 117 membrane was immersed in 5% H 2 O 2 solution at 80°C for 1 h, then in 0.5 M H 2 SO 4 solution at 80°C for an additional hour, and finally washed with deionized water. The electrode preparation and electrochemical measurements are done using the previously reported procedure 55 . An Ag/AgCl electrode (3.5 M KCl) and a Pt foil were used as the reference and counter electrodes, respectively. The potentials were measured against the Ag/AgCl electrode and converted to RHE according to E (vs. RHE) = E (vs. Ag/AgCl) + 0.198 V + 0.059 × pH. The working electrodes were prepared as follows: 5 mg of carbon-supported catalyst, 1 ml isopropanol, and 20 µl of Nafion solution (5 wt%, Sigma-Aldrich) were mixed and ultrasonicated for more than 30 min to generate a homogeneous ink. Then, 200 µl of the catalyst ink was deposited onto a 1 cm 2 carbon fiber paper by drop-casting, resulting in a metal loading of~0.2 mg cm −2 . Before the electrochemical measurement, the 1 M KOH + 1 M KNO 3 electrolyte was purged with Ar for at least 30 min. The LSV curves were obtained by scanning the potential from −0.5 to −1.7 V vs. Ag/AgCl at a rate of 20 mV s −1 . CA measurements were conducted at various potentials in 1 M KOH + 1 M KNO 3 solution with an Ar flow rate of 20 sccm in the cathodic compartment. CV measurements at various scan rates were performed to estimate the C dl of the catalysts in a static solution by scanning the potential between −0.1 and 0 V vs. Ag/AgCl. EIS measurements were carried out in the potentiostatic mode at −1.2 V vs. Ag/AgCl, applying a 5 mV AC dither and scanning from 100 kHz to 100 MHz.
Quantification of products. The gas product was quantified by gas chromatography (Agilent 7890B). The quantity of ammonia produced was measured using a colorimetric method with Nessler's reagent. All test solutions were incubated under dark conditions at room temperature for 20 min before UV-vis tests. The absorbance at 420 nm for each solution was measured with a UV-vis spectrophotometer (Agilent 3500). A series of reference solutions with suitable NH 4 Cl concentrations was created to plot a calibration curve. Electrolytes after catalysis were diluted to ensure the ammonia concentrations in the test solutions were in the linear range of Nessler's method. The concentrations of ammonia in the electrolytes were obtained with this as-obtained calibration curve. Ion chromatography was also used to detect the liquid products. The liquid products from the isotopic experiment were analyzed by 1 H NMR using dimethyl sulfoxide (DMSO) (20%) and 1 mM maleic acid as the internal standard. The pH of the after-reaction electrolyte was adjusted to 2 using 1 M HCl. Faradaic efficiency was calculated according to the following equation: Yield rate was calculated according to the following equation: where n is the number of electron transfers towards the formation of 1 mol product; C is the concentration of product (M); V is the volume of catholyte (mL); F is the Faraday constant (96,485 C • mol −1 ); i is the reduction current; t is the total reaction time, and m is the catalyst mass.
DFT calculations. DFT calculations were performed using VASP 56,57 and Quantum ESPRESSO (QE) 58 at the GGA level using the RPBE 59 functional. QE calculations were used to get the *N and *NO 3 adsorption energies and the atom and molecular projected density of states used to optimize the Bayeschem models. For QE calculations, the core electrons were treated using ultrasoft pseudopotentials with kinetic energy cutoffs of 500 eV for wavefunction and 5000 eV for charge density. In order to speed up calculations, partial occupancies were set using the Fermi-Dirac smearing with a smearing parameter of 0.1 eV. For all the other electronic structure calculations, VASP with projector augmented wave pseudopotentials was used. A planewave energy cutoff of 450 eV was used and the same kpoint setting as QE were used for VASP. The Methfessel-Paxton smearing scheme was used with a smearing parameter of 0.2 eV for adsorbate systems and 0.001 eV for molecules. Electronic energies are extrapolated to k B T = 0 eV. For both VASP and QE the same model systems were used. The (100) surfaces were modeled using a 3 × 3 × 6 supercell with the bottom 4 layers fixed, and the (111) surfaces were modeled using a 4 × 4 × 4 supercell with the bottom two layers fixed. The Brillouin zones for (111) and (100) surfaces were sampled using the Monkhorst-Pack meshes of 3 × 3 × 1 and 4 × 4 × 1, respectively. A vacuum of 15 Å was used between two periodic images. Geometries were optimized until the maximum forces was less than 0.05 eV/Å. Calculations for molecules or atoms were done in a 15 × 15 × 15 Å box using the Gamma point. Atom and molecular projected density of states 60 were calculated in QE on a denser k-point sampling of 12 × 12 × 1 with an energy spacing of 0.01 eV.
To convert the DFT electronic energies to free energies, we added the zeropoint energy (ZPE), heat capacity, and entropic contributions. ZPE and entropic corrections were calculated within the harmonic oscillator approximation. For adsorbates, all degrees of freedom are treated as vibrational. Using the Schumann et al. 61 approach, frequencies less than 100 cm −1 are treated as pseudotranslational/rotational and thus replaced by 12 cm −1 . For gas-phase species, translational and rotational contributions to the internal energy and entropy are considered using statistical thermodynamics. All the corrections are calculated at 298 K and tabulated in Table S1. The free energy of liquid water is calculated using the gas phase at its vapor pressure at 298 K (0.035 bar). The free energy of an adsorbed system *A at the gas-solid interface is given by the following equation: At the solid-electrolyte interface under constant potential conditions, the Gibbs free energy of an adsorbed system *A is then calculated as where μ e is the chemical potential of the electron (same as the Fermi level since the vacuum potential is 0). N e is the number of electrons added (positive in sign) or removed (negative in sign) relative to the total number of electrons in the chargeneutral system. In this case, E GCÀDFT ÃA is the electronic energy calculated from the GC-DFT. Grand-canonical DFT calculations (GC-DFT) were performed in order to simulate the solid-electrolyte interfaces with the applied bias and solvation. In GC-DFT calculations, the applied potential bias is controlled by changing the number of electrons in the simulation cell, which modifies the Fermi level relative to the vacuum potential and thus the work function of the system (ϵ F = −ϕ). The equation below shows how the work function (ϕ) is connected to the applied potential at the SHE (standard hydrogen electrode) scale.
Therefore, systems can be simulated at a specific potential by optimizing the number of electrons to reach the target work function. Solvation was included through the VASPsol 62,63 with a continuum dielectric description of electrolytes. The surface tension parameter (0), the dielectric constant (78.4), and the Debye screening length (3 Å) were set in all GC-DFT calculations 64,65 . The electrode potential correction to the free formation energy of a surface intermediate was obtained from GC-DFT calculations on Cu(100) at 0 V vs RHE. This correction is assumed to be constant, which is a reasonable approximation across metal surfaces at this potential (Supplementary Table 2). Examples of calculating adsorbate free formation energies are shown in the Supplementary discussion and thermodynamic cycle in Supplementary Fig. 11 for nitrate adsorption. We include the most stable adsorbate configurations for pure metals in a Github repository (https://github.com/hlxin/nitraterr).
Bayesian theory of chemisorption (Bayeschem). Within the d-band theory of chemisorption, the adsorption energy (ΔE A ) of a species at metal surfaces can be calculated as: where ΔE 0 is the energy change due to the interaction between the adsorbate and the free-electron-like sp-states of the transition metals, while the ΔE d is due to the interaction between the adsorbate and the localized d-states. It is assumed that the sp-band contribution is constant for a given surface facet due to similar sp-states of the transition metals. Whereas the variations in binding energies are explained due to differences in the d-contribution 66 . The d-contribution can be further decomposed into repulsive orbital orthogonalization and attractive orbital hybridization.
The d-hybridization contribution can be calculated with the Newns-Anderson model, which depends on several parameters (α, β, Δ 0 , ϵ a ), as shown below for a simple case with one valence state of the adsorbate.
where the chemisorption function Δ(ϵ) depends on the surface density of states (ρ d ): The repulsive contribution from each adsorbate frontier orbital is derived by treating the adsorbate-substrate interaction as a two-level system.
in which f and e n a are the filling of the metal d-states and the adsorbate resonance state. A Bayesian learning approach was used to optimize the model parameters mentioned above by learning from ab initio adsorption properties. Further details on this approach can be seen in ref. 15 . Bayesian models were optimized for *NO 3 and *N on the fcc(100) and fcc (111)

Data availability
The data that support the plots within this paper and other findings of this study are available from the corresponding author upon reasonable request.